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Abstract 

Ferroelectric  materials,  such  as  PZT,  PLZT,  PMN  and  BaTiC>3,  provide  unique  actuator  and 
sensor  capabilities  for  applications  including  nanopositioning,  high  speed  valves  and  fuel  injectors, 
camera  focusing  and  shutter  mechanisms,  ultrasonic  devices  for  biomedical  imaging  and  treatment, 
and  energy  harvesting  devices.  However,  to  achieve  the  full  potential  of  the  materials,  it  is  necessary 
to  develop  and  employ  models  that  quantify  the  creep,  rate-dependent  hysteresis,  and  constitutive 
nonlinearities  that  are  intrinsic  to  the  materials  due  to  their  domain  structure.  The  success  of  mod¬ 
els  requires  that  they  be  highly  efficient  to  implement  since  real-time  applications  can  require  kHz 
to  MHz  rates.  The  calibration  of  models  for  specific  materials,  devices,  and  applications,  requires 
efficient  and  robust  parameter  estimation  algorithms.  Finally,  control  designs  can  be  facilitated 
by  models  that  admit  efficient  and  robust  approximate  inversion.  The  homogenized  energy  model 
(HEM)  is  a  multiscale,  micromechanical  framework  that  quantifies  a  range  of  hysteretic  phenomena 
intrinsic  to  ferroelectric,  ferromagnetic  and  ferroelastic  materials.  In  this  paper,  we  present  highly 
efficient  implementation  and  parameter  estimation  algorithms  for  the  ferroelectric  model.  This  in¬ 
cludes  techniques  to  construct  analytic  Jacobians  and  data-driven  algorithms  to  determine  initial 
parameter  estimates  to  facilitate  subsequent  optimization.  The  efficiency  of  these  algorithms  fa¬ 
cilitates  material  and  device  characterization  and  provides  the  basis  for  constructing  efficient  and 
robust  inverse  algorithms  for  model-based  control  design.  The  model  implementation,  calibration, 
and  validation  are  illustrated  using  rate-dependent  PZT  data  and  single  crystal  BaTiC>3  data. 
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Nomenclature 


d  Piezoelectric  constant  (m/V  =  C/N) 

h  Piezoelectric  constant  (V/m  =  N/C) 

E  Electric  field  (V / nr) 

g  Gibbs  energy  for  dipoles  (CV) 

Ga  Gibbs  energy  density  of  a- variant  (CV/nr3) 

Gea  Electric  Gibbs  energy  density  of  a- variant  (CV/nr3) 

P  Polarization  (C/nr2) 

P  Polarization  kernel  (C/nr2) 

P°  Polarization  of  cr- variant  (C/nr2) 

Pr  Remanent  polarization  of  a- variant  (C/nr2) 

Minimum  polarization  of  a- variant  (C/nr2) 

Elastic  compliance  of  a-variant  at  constant  field  (nr2/N) 

T  Temperature  (K) 

x_|_,x_,xgo  Fraction  of  positively,  negatively  and  90°  oriented  dipoles  (Unitless) 

Elastic  stiffness  of  a-variant  at  constant  polarization  (N/nr2) 
e  Permittivity  (F/nr  =  C/Vnr) 

e  Strain  (Unitless) 

e  Strain  kernel 

ea  Strain  of  a-variant 

Renranence  strain  of  a-variant 
Minimum  strain  of  a-variant 

if  Inverse  susceptibility  at  constant  strain  (nr/F  =  Vnr/C) 

a  Stress  (N/nr2) 

T90  Relaxation  time  for  90°  switching  (s) 

Tiso  Relaxation  time  for  180°  switching  (s) 

Xe  Electric  susceptibility  (Unitless) 

Xa  Ferroelectric  susceptibility  of  a-variant  at  constant  stress  (F/nr  =  C/Vnr) 

ipa  Helmholtz  energy  density  of  a-variant  (CV/nr3) 

1  Introduction 

Ferroelectric  materials,  such  as  lead  zirconate  titanate  (PZT),  lanthanunr-doped  lead  zirconate  ti- 
tanate  (PLZT),  lead  manganese  niobate  (PMN),  polyvinylidene  flouride  (PVDF),  and  barium  ti¬ 
tanate  (BaTiOs),  are  being  widely  considered  as  actuators,  sensors,  and  structural  units  for  a  range 
of  aerospace,  aeronautic,  automotive,  industrial  and  biomedical  applications.  Their  advantages  derive 
from  the  intrinsic  properties  of  the  materials.  The  complementary  direct  and  converse  piezoelectric 
effects  imbue  the  materials  with  multiple  design  properties  including  sensor,  actuator,  self-monitoring 
and  nondestructive  evaluation,  and  energy  harvesting  capabilities.  Their  functionality  is  augmented 
by  their  solid  state  nature  which  promotes  miniaturization  and  simplified  design  and  reduces  power 
requirements  and  heat  generated  by  the  units.  Furthermore,  ferroelectric  actuators  can  provide 
nanometer  positioning  resolution  and  can  operate  at  frequencies  ranging  from  DC  to  MHz.  As  de¬ 
tailed  in  the  companion  paper  [19],  ferroelectric  transducers  are  being  considered,  or  are  already 
being  employed,  for  a  large  number  of  applications  including  high  speed  valves  for  fuel  injectors,  fer¬ 
roelectric  memory  technologies  (e.g.,  FeRAM),  nanopositioning  units  (e.g.,  AFM  and  STM  stages), 
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high  speed  camera  shutters  and  autofocusing  units,  ultrasonic  transducers  for  biomedical  imaging 
and  treatment,  micro  and  pico  air  vehicle  design,  and  energy  harvesting  devices. 

However,  the  ferroelectric  properties  that  imbue  the  materials  with  unique  transducer  capabil¬ 
ities  also  produce  rate-dependent  hysteresis,  creep,  and  constitutive  nonlinearities  as  illustrated  in 
Figure  1.  The  field-polarization  and  field-strain  data  from  [26]  illustrates  that  rate-dependent  effects 
are  significant  at  frequencies  as  low  as  1  Hz.  PZT  data  from  [24]  illustrates  that  for  fixed  field  inputs, 
both  the  strain  and  polarization  exhibit  significant  creep  on  timescales  of  1  to  20  seconds.  Finally, 
the  MFC  data  from  [8]  illustrates  nested  minor  loop  behavior  typical  of  moderate  drive  regimes. 
Whereas  the  full  switching  behavior  shown  in  Figure  l(a)-(e)  will  typically  not  be  encountered  in 
applications,  general  models  must  account  for  the  full  range  of  rate-dependent  and  creep  behavior 
to  provide  comprehensive  device  characterization. 

To  be  optimally  utilized  in  applications,  constitutive  models  must  satisfy  the  following  criteria: 

(i)  they  must  adequately  quantify  the  range  of  rate-dependent  behavior  exhibited  by  the  materials, 

(ii)  they  must  be  in  a  form  that  can  be  readily  extended  to  provided  distributed  models,  for  complex 
structures  and  devices,  that  are  amenable  to  finite  element  implementation,  (iii)  they  must  be  efficient 
to  implement,  (iv)  they  must  be  readily  calibrated  for  a  specific  material  or  device,  and  (v)  models 
employed  in  control  designs  can  prove  advantageous  if  they  admit  the  construction  of  approximate 
inverse  relations. 

It  is  shown  in  [18-21]  that  the  homogenized  energy  model  (HEM)  satisfies  (i)  and  (ii).  Fur¬ 
thermore,  the  construction  of  approximate  inverses  for  previous  formulations  of  the  framework  are 
reported  in  [4,18]  and  inversion  techniques  for  the  present  formulation  [19]  are  addressed  in  [15].  In 
this  paper,  we  address  (iii)  and  (iv)  by  presenting  highly  efficient  implementation  techniques  and 
data-driven  parameter  estimation  algorithms  for  the  homogenized  energy  model.  In  combination,  this 
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Figure  1:  (a)  Field-polarization,  (b)  field-strain  and  (c)  time-strain  PZT  data  from  [24],  Rate- 
dependent  (d)  polarization  and  (e)  strain  PZT  data  from  [26].  (f)  Field-strain  MFC  data  from  [8]. 
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provides  the  framework  with  significant  flexibility  for  use  in  simulation  packages,  design  algorithms, 
and  control  designs  for  systems  utilizing  ferroelectric  transducers. 

To  provide  context  for  the  framework,  we  first  summarize  the  mechanisms  that  produce  the  hys- 
teretic  behavior  shown  in  Figure  1.  The  creep  and  rate-dependencies  are  due  to  the  fact  that  kinetics 
associated  with  dipole  switching  typically  differ  from  mechanical  or  electrical  loading  rates.  Further¬ 
more,  it  is  discussed  in  [19]  that  the  180°  switching  associated  with  large  polarization  changes  often 
occurs  at  a  different  rate  than  90°  switching  associated  with  large  strain  changes.  In  combination, 
this  establishes  that  multiple  time  scales  must  be  incorporated  in  dynamic  models. 

As  illustrated  for  a  PZT-based  macro-fiber  composite  (MFC)  in  Figure  2,  hysteresis  also  involves 
multiple  spatial  scales.  Within  a  domain,  strains  or  changes  in  polarization  are  due  to  stress-induced 
material  deformations  or  field-induced  ion  movement  and  the  behavior  is  often  reversible  and  linear. 
Field  or  stress-induced  dipole  switching  at  the  grain  level  produces  irreversible  hysteresis  in  both  the 
field-polarization  and  field-strain  relations.  For  single  crystal  materials  comprised  of  a  single  grain, 
the  switching  is  typically  rapid  thus  producing  sharp  hysteresis  and  butterfly  loops  in  quasistatic 
operating  regimes.  For  polycrystalline  materials  with  distributed  interaction  and  coercive  fields, 
hysteresis  and  butterfly  loops  are  smoothed  due  to  nonuniform  grain  contributions.  The  behavior 
of  hysteresis  loops  is  further  modified  when  hysteretic  actuator  or  sensor  materials  are  employed  on 
distributed  structures. 


Imerdigitated 

Electrodes 


Material  Level 

P(E,  a)  =  d(E,  a)a  +  X°E  +  Pirr(E,  a) 
e{E,  a)  =  sEa  +  d(E,  a)E  +  £irr(E,  a) 


Grain  Level 

P(E,a)  =  xaE  +  £(d0< 7  +  P^)xa{E,a) 

a 

HE,  a)  =  sEa  +  ^2,(daE  +  eaR)xa{E ,  a) 


Figure  2:  Multiscale  behavior  of  a  PZT-based  MFC  transducer  at  the  application,  device,  material 
and  unit  cell  levels. 
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As  noted  in  Figure  2,  material-level  constitutive  relations  for  the  polarization  P  and  strain  e  can 
generically  be  expressed  as 


P(E,  a)  =  d(E,  a) a  +  x°  E  +  Pirr(E,  a) 
s(E,  a)  =  sEa  +  d(E,  a)E  +  £irr{E ,  a) 

where  E,  cr  are  input  fields  and  stresses,  xa  is  the  ferroelectric  susceptibility  at  constant  stress,  sE 
is  the  elastic  compliance  at  constant  field,  and  d  is  a  piezoelectric  coupling  coefficient.  We  note  that 
d(E,a),  Pirr(E,a)  and  £irr(E,a)  incorporate  the  nonlinear  and  irreversible  history-dependence  due 
to  dipole  switching. 

Various  modeling  hierarchies  can  be  defined  by  the  manner  in  which  d,  Pirr  and  e,rr  are  con¬ 
structed.  Micromechanical,  or  microscopically-motivated  models  are  based  on  an  energy  description 
of  the  material  at  the  domain  or  grain  level  in  combination  with  various  homogenization  techniques  to 
provide  expressions  for  the  nonlinear,  effective  components  d,Pirr  and  £jrr.  Phenomenological  mod¬ 
els  circumvent  the  difficulties  associated  with  quantifying  complex,  or  poorly  understood,  micro-scale 
physics  by  constructing  relations  for  d,  Pirr  and  £irr  based  on  macroscale  observations  or  experimen¬ 
tal  measurements,  often  guided  by  thermodynamic  constraints.  Details  regarding  various  models  for 
ferroelectric  materials  can  be  found  in  [2, 9, 13, 16, 18, 19]. 

The  homogenized  energy  model  is  a  multiscale,  nricronrechanical  approach  that  begins  with  energy 
analysis  at  the  domain  level  to  construct  the  linear  relations  shown  in  Figure  2  for  Pa  and  £a  where 
a  designates  the  dipole  variant;  e.g.,  a  =  ±180,90  for  tetragonal  materials.  Switching  processes 
due  to  domain  wall  nucleation  and  movement  are  incorporated  by  tracking  the  evolution  of  dipole 
fractions  xQ  which  serve  as  internal  variables.  Differential  equations  quantifying  the  dynamics  of  xa 
are  driven  by  likelihood  rates  constructed  using  Boltzmann  theory  to  quantify  the  scaled  probability 
of  transitioning  between  stable  equilibria  associated  with  dipole  variants.  This  incorporates  the 
rate-dependence  and  multiple-time  scales  exhibited  by  the  data  in  Figure  1  and  yields  grain-level  or 
single  crystal  kernels  P  and  e.  For  poly  crystalline  materials,  material-level  models  are  constructed 
by  assuming  that  properties  such  as  coercive  fields,  critical  driving  forces,  and  interaction  fields  are 
manifestations  of  underlying  densities  rather  than  constants. 

The  three  primary  issues  that  must  be  addressed  to  construct  implementation  algorithms  are:  (i) 
efficient  approximation  of  the  integrals,  posed  on  infinite  and  semi-infinite  domains,  arising  in  the 
homogenization  step  used  to  construct  the  irreversible  components  d(E,  a),  Pirr(E1  a)  and  £irr(E ,  a), 
(ii)  efficient  approximation  of  the  evolution  equations  for  xa,  and  (iii)  efficient  evaluation  of  the  like¬ 
lihoods  which  incorporate  rate-dependent  effects  and  drive  the  evolution  equations.  The  first  issue 
is  addressed  by  employing  density  representations  that  admit  efficient  numerical  approximation  on 
finite  domains  whereas  implicit  Euler  discretizations  readily  address  the  second  issue.  The  repeated 
evaluation  of  likelihoods  is  avoided  by  constructing  arrays  that  permit  high-speed  access  during  model 
implementation.  In  combination,  the  evaluation  of  d(E,a),Pirr(E,a )  and  £irr(E,a)  are  reduced  to 
operations  solely  involving  componentwise  matrix  multiplication  and  summation.  Hence  the  oper¬ 
ations  are  highly  efficient  and  inherently  parallelizable  which  facilitates  subsequent  implementation 
using  devices  such  as  field  programmable  gate  arrays  (FPGA)  [3] . 

We  note  that  approximate  model  inversion  techniques  rely  on  forward  solutions  of  the  model. 
Hence  the  efficiency  of  forward  algorithms  also  directly  contributes  to  the  speed  of  inverse  algorithms. 

There  are  two  primary  issues  that  must  be  addressed  when  employing  gradient-based  optimization 
algorithms  for  the  parameter  estimation  required  for  model  calibration:  (i)  determination  of  accurate 
initial  parameter  estimates  and  (ii)  efficient  and  accurate  construction  of  gradients  or  Jacobians.  Due 
to  its  physical  nature,  parameters  in  the  homogenized  energy  model  can  be  directly  correlated  with 
properties  of  measured  data.  We  utilize  that  here  to  construct  data-driven  algorithms  to  determine 
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initial  parameter  estimates.  Secondly,  we  construct  analytic  Jacobian  relations  which  significantly 
improves  both  the  speed  and  accuracy  of  optimization  routines.  In  combination,  this  helps  mitigate 
some  of  the  difficulties  associated  with  the  shallow  slopes  and  nonconvexity  (multiple  local  minima) 
inherent  to  functionals  employed  in  parameter  estimation  problems  of  this  type. 

We  note  that  various  aspects  of  the  implementation  and  data-driven  parameter  estimation  algo¬ 
rithms  can  be  optimized  to  provide  improved  performance  for  specific  data  sets.  Instead,  our  goal 
was  to  provide  algorithms  that  balance  simplicity,  robustness,  accuracy,  and  optimality  for  a  wide 
range  of  materials,  applications  and  operating  regimes.  Points  were  algorithms  can  be  modified  or 
optimized  are  noted  in  the  discussion. 

Following  a  brief  summary  of  the  model  in  Section  2,  highly  efficient  implementation  algorithms 
are  presented  in  Section  3.  The  parameter  estimation  problem  is  discussed  in  Section  4  where  we 
provide  algorithms  to  construct  analytic  Jacobians  and  data-driven  techniques  to  determine  initial 
parameter  estimates.  Model  calibration  and  validation  are  illustrated  in  Section  5  using  PZT  data 
from  [24]  and  single  crystal  BaTiC>3  data  from  [5].  Additional  examples  illustrating  the  capability  of 
the  model  to  characterize  creep  and  rate-dependent  data,  such  as  the  shown  in  Figure  1,  are  provided 
in  the  companion  paper  [19]  which  details  the  model  development. 


2  Polarization  and  Strain  Models 


As  detailed  in  [19],  consideration  of  180°  dipole  switching  yields  a  macroscopic  model  quantifying 
the  nonlinear  and  hysteretic  field-polarization  map  for  ferroelectric  materials  whereas  the  additional 
incorporation  of  90°  switching  yields  a  homogenized  model  quantifying  both  the  polarization  and 
strains  due  to  input  fields  and  stresses.  We  summarize  the  latter  model  and  indicate  simplifications 
that  occur  if  considering  only  the  polarization  in  the  absence  of  applied  stresses. 


Polarization  and  Strain  Model 

At  the  lattice  level,  we  consider  the  Helmholtz  and  Gibbs  energy  densities 

1  „  _ o  1 


MP,e)  =  2 ~  PrY  +  2Va  >  -  £rY  +  ha(P  ~  Pg)(e  -  s%) 


and 


Ga{E,  a;  P,  e)  =  i/j a(P ,  e)  —  EP  —  ae 

where  we  indicate  ±180°  and  90°  orientations  by  a  =  ±,  90.  As  summarized  in  the  nomenclature 
table  at  the  beginning  of  the  paper,  Pr,£%  Yp  and  ha  denote  the  remanence  polarization,  rerna- 
nence  strain,  inverse  susceptibility  at  constant  strain,  elastic  stiffness  at  constant  polarization,  and 
piezoelectric  constant. 

For  a  fixed  applied  field  and  stress,  the  conditions  =  0  and  =  0  yield  the  relations 


PY  =  Pr  +  XaaE  +  dQa 
£m  =  £r  +  daE  +  s^a 


(1) 


where 


'X.a 


Yp 

a 


i  da  — 


hn 


sE  = 


To 


Ypne  -  h2  ’  “  h 2  -  Ypn£  ’  “  Ypn£  -  h2 

1  a.  'la.  ,La  ILa  1  a  'la  1  a  'la  ILa 

are  the  ferroelectric  susceptibility  at  constant  stress,  the  piezoelectric  constant,  and  elastic  compliance 
at  constant  field.  The  minimum  of  the  Gibbs  energy  in  each  a-well  can  then  be  expressed  as 

Gam(E,  a)  =  -\x*E2  -  ^ a2  -  daEo  -  EP%  -  ae% 
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The  dipole  fractions  x+,x-  and  X99  associated  with  positively,  negatively,  and  90°  degree  dipoles 
evolve  according  to  the  differential  equations 

x-  =  (P—90  +  P-+)x-  +  P90-X90  +  V+-X+ 

±99  =  P-99X-  -  (p90-  +  P90+)®90  +  P+ 90X+  (2) 

X+  =  P-+X-  +P90+X99  -  (p+90  +p+-)x+ 


where 


Pap{E,cr) 


1  c-AG“a(E,o)V/kT 
Taf3 


(3) 


quantifies  the  likelihood  of  transitioning  from  an  a- well  to  a  /3-well.  As  noted  in  Table  1,  (2)  can  be 
simplified  using  the  relation  x+  +  X-  +  xgo  =  1 .  The  activation  energy  is  specified  by  the  relation 


A  Gaa/3(E,a;Fc) 


AG0(1  -  Fap(E,  a)/Fc)2  ,  Fap(E,  a)  <  Fc 
0  ,  Fap(E,  a)  >  Fc. 


Here  Fap(E,  a)  =  Gam(E,  a)  —  Gpm(E,  a)  is  the  thermodynamic  driving  force.  The  specific  form  of 
Fap,  based  on  the  physical  assumption  that 


X+  =  X-  =  X90  =  X 

S+  =  S-  =  S90  =  SE 

Fr  =  0)  Pr  =  -Pr  >  eR  =  eR, 

dgo  =  0  ,d-  =  -d+ 


(4) 


790—  —  r-90  —  T90+  —  7+90  —  Tgo  ,  T+_  —  T_+  —  Ti80, 


is  given  in  Table  1.  Note  that 


AG0  = 


-  F 

4  r° 

ip 

16Xc 


180°  Switching 
90°  Switching 


is  the  energy  barrier  at  zero  driving  force. 

The  polarization  and  strain  kernels  are  given  by 


p=  E  x«p“ 

a=±,90 


£  = 


E 


Xa£ 


a 

m* 


a=±,90 


Based  on  (1)  and  (4),  these  relations  can  be  expressed  as 

P{E,  a)  =  d(E ,  a)a  +  xaE  +  Pirr{E,  u) 
e(E,  a)  =  sEcr  +  d(E ,  a)E  +  £irr(E,  a) 


where 

d(E,a)=  ^2  daxa(E,cr ) 

a=±,90 

P irr(E,  (7 )  =  ^  P%Xa(E,C r) 

a=±,90 

£irr(E,a)=  ^2  £RXa(E,a). 
a=±,90 
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In  the  final  step  in  the  development  of  [19],  macroscopic  models 

/*oo  /»oo 

P(E(t),a(t)-x°+)=  /  PiEW  +  E^aWiFcWErivciFjdEjdFc 

J  0  J — oo 

/*oo  /»oo 

£(E(t),a(t)-x°+)=  /  E(E(t)  +  EI,a(t)-,Fc)vI(EI)i'c(Fc)dEIdFc 

J  0  J — oo 

are  constructed  by  considering  interaction  fields  £/  and  thermodynamic  driving  forces  Fc  to  be 
manifestations  of  underlying  densities  uj  and  vc.  The  final  models  can  thus  be  expressed  as 

P(E ,  a)  =  d(E,  a) a  +  P  +  Pirr{E,  a) 
e(E,  a)  =  sEa  +  d(E,  a)E  +  £irr(E,  a) 


where 

/*oo  /'OO 

d(E,a)  =  /  /  d(Ee-,Fc)uI(EI)iyc(Fc)dEIdFc 

J  0  J —oo 

r  oo  /*oo 

Pirr(E,  a)  =  I  /  Pirr(Ee;Fc)vI(EI)i'c(Fc)dEIdFc 
J  0  J — oo 

/»oo  /*oo 

—  /  I  £irr(Ee]  FC^V] (Ei^lSc^Fo^dEjdFc. 

J  0  J  —oo 

Here  Ee(t)  =  E(t )  +  £/  is  the  effective  field. 

Details  regarding  various  choices  for  the  densities  are  provided  in  [19].  For  the  implementation 
and  parameter  estimation  algorithms  detailed  in  this  paper,  we  employ  the  representations 


Ay, 


vc{Fc)  =  akMFc )  ,  &(FC)  =  - - - .-KAc)-in(Fc)]V2(^) 


i 


k - /Cq; 


a^Fcy/2iT 


MEi)  =  vTt  X!  Pk<Pk{Ei)  ,  <^(£7)  = 


2  fry/ „k\2  k  _  u 

i  ° c  —  2  Uc 


_fc  nA; 

,  aI  =  2  <n 


(5) 


where  the  preceding  sums  ensure  integration  to  unity.  During  model  calibration,  the  parameters 
{ak,Pk}  are  determined  through  a  least  squares  fit  to  the  data. 

The  complete  polarization  and  strain  model  is  summarized  in  Table  1. 


Parameters  for  the  Polarization-Strain  Model 

To  construct  the  density  basis  functions  ( pk{Ei )  and  <f>k(Fc),  it  is  necessary  to  determine  values 
for  the  driving  force  mean  [ic  =  ln(£c)  and  variance  cr2  as  well  as  the  interaction  field  variance  aj. 
We  denote  this  set  of  parameters  by 

P=\Mc,o%,oj\-  (6) 

We  provide  an  algorithm  in  Section  4.2  to  estimate  values  for  these  parameters  based  on  measured 
attributes  of  E-P  data.  We  note  that  these  parameters  are  not  updated  or  optimized  through  a  least 
squares  fit  to  data. 

Based  on  the  assumption  (4),  the  remaining  parameters  are  denoted  by 

P=  [PR,£R^R,X(7,d+,sE,'y,T90,T1so,ak,f3k\.  (7) 

This  comprises  the  set  that  is  optimized  during  model  calibration. 
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Homogenized  Energy  Model: 


where 


P{E,  a)  =  d(E,  a)  cr  ±  \a E  +  Pirr(E,  a) 
e(E,  a)  =  -sEa  ±  d(E,  a)E  ±  £irr(E,  a) 

poo  poo 

d{E,  a )  =  /  /  d(E  +  Ej\  Fc)yj{Ej)uc(Fc)dEjdFc 

J  0  J  —  oo 

poo  poo 

Pirr(E,a)  =  /  /  Pirr(E  +  EI;Fc)vI(EI)isc(Fc)dEIdFc 

J  0  J —oo 

poo  poo 

£irr{E,a)  =  /  /  £irr(E  +  EI-,Fc)i^I(EI)iyc(Fc)dEIdFc 

J  0  J —oo 


Kernels: 


d(E,  a;  Fc)  =  y  daxa(E;Fc ) 

a=d=,90 

Pirr(E,  a;  Fc)  =  J2  PrMF;Fc) 
a=±,90 

£irr{E,a;Fc)  =  y  £%xa(E;Fc) 

a=±,90 

Evolution  Equations:  (t;  E.  <7,  Fc) 

x-  =  ~{P-90  +  P-  +  +P9 0-)X-  ±  (p+-  ~P90-)x+  +P99- 
x+  =  0-+  ~P90+)X-  -  {p+90  +  P+-  +P90+)x+  +P90+ 
X90  =  1  —  X+  —  X- 

Likelihoods:  a,  (3  =  ±90 


Densities: 


Ka 


1  *  1 

Vc(Fc)  —  ^  '  Qk4>k{.Fc)  ,  (pk^Eoj 

z2eat 


3-[ln(Fc)-ln(Fc)]2/2(^)2 


1 


k=ka 


vi{Ei)  =  Y  f3kVk(Ei)  ,  (fk{Ej)  = 


a^Fcy/2 7r 

e-Ef/2(^)2 


(8) 


(9) 


(10) 


(11) 


Xa0 

7  =  V/feT 

(12) 

Activation  Energy:  a,  (3  =  ±90 

,  Fap(E,  a)  <  Fc 
,  Fap(E,  a)  >  Fc 

F90+(E,  a)  =  d+Eo  ±  P+E  -  (4°  -  £+)a 

1  E+90  =  -F90+ 

(13) 

E9o-{E,a)  =  -d+Ea  -  P+E  -  (4°  -  e+)a 

,  E-  90  =  —E90— 

F _| _ (E,  a)  =  —2d+Ea  —  2 EP^ 

1 

II 

+ 

(14) 


Table  1:  Components  of  the  polarization  and  strain  model  based  on  90°  dipole  switching. 


Polarization  Model 

If  solely  quantifying  the  polarization  in  the  absence  of  applied  stresses,  one  can  consider  energy 
landscapes  and  dipole  fractions  associated  only  with  positively  and  negatively  oriented  dipoles  which 
we  designate  by  a  =  ±.  It  is  shown  in  [19]  that  a  complete  characterization  of  the  Helmholtz  and 

Gibbs  energy  densities  yields  likelihoods  p and  p. \ formulated  in  terms  of  complementary  error 

functions  whereas  formulation  in  terms  of  the  activation  energy  yields  likelihood  relations  analogous 
to  (3).  We  focus  on  this  latter  case  and  note  that  analogous  implementation  and  parameter  estimation 
algorithms  can  be  constructed  for  the  error  function  formulation. 

The  polarization  model  based  on  180°  dipole  switching  is  summarized  in  Table  2. 

Polarization  Model  Parameters 

In  this  case,  we  employ  a  coercive  field  density  vc(Ec)  having  the  same  form  as  the  driving  force 
density  vc{F^).  Hence  it  is  again  necessary  to  prescribe  data-driven  algorithms  to  determine  values 
of  the  coefficients  p  in  (6).  The  parameters  to  be  estimated  through  a  least  squares  fit  to  data  are 

P=[v,pR,rY,T,ak,/3k\.  (15) 

We  discuss  associated  optimization  algorithms  in  Section  4  and  data-driven  techniques  to  obtain 
initial  parameter  estimates  in  Section  4.2. 


3  Implementation  Algorithms 

To  facilitate  implementation  of  the  models  summarized  in  Tables  1  and  2,  or  their  inverses,  three 
issues  must  be  addressed:  (i)  approximation  of  the  infinite-domain  integrals  in  a  manner  that  retains 
its  accuracy  for  the  basis  expansions  (5),  (ii)  evaluations  of  the  likelihoods  pan  in  an  efficient  manner, 
and  (iii)  efficient  solution  of  the  differential  equations  (11)  and  (19). 


3.1  Quadrature  Techniques 

To  implement  the  model,  the  integrals  defining  the  terms  d,  Pirr  and  £irr  in  (9)  and  (17)  must  be 
approximated  in  a  manner  that  is  both  efficient  and  accurate.  We  first  note  that  by  employing  the 
relations  (10),  these  terms  can  be  expressed  as 


d(E,a)  = 

Pirr(E ,  (j) 
&irr(E ,  a) 


£  4 

q=±,90 


=  £  n 

a=±,  90 


xa{Ee ;  Fc)fI{EI)uc{Fc)dEIdFc 


poo 

I  xa(Ee-,Fc)i/i(EI)i/c(Fc)dEIdFc 

—  OO 

•oo 


=  y}2£R  /  xa(Ee]  Fc)iyI(EI)uc(Fc)dEIdFc. 

n — t-  an  do  J —oo 


Secondly,  it  is  noted  in  [19]  that  the  densities  satisfy  exponential  decay  constraints.  Hence  approxi¬ 
mation  algorithms  can  be  defined  on  finite  domains  rather  than  necessitating  quadrature  techniques 
for  infinite  and  semi-infinite  domains.  We  illustrate  here  a  trapezoidal  formula  and  refer  the  reader 
to  [18]  for  details  regarding  composite  Gaussian  quadrature  techniques.  We  note  that  truncation 
of  the  domains  and  use  of  composite  techniques  is  more  efficient  to  implement  than  Gauss-Hermite 
and  Gauss-Laguerre  algorithms,  that  are  designed  specifically  for  infinite  and  semi-infinite  intervals, 
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Homogenized  Energy  Model: 


P(E)  — - b  Pirr 

n 


where 


Kernel: 


Pit 


/*oo  /*oo 

■r(E)  =  /  /  Pi.rr{E  +  Ej;  Fc)iyj(Ej)uc(Fc)dEjdFc 

J  0  J  — oo 

Pirr(E]  Fc)  =  PftXa(E',  Fc)  =  —  Hr  +  2Prx_|- 


a=±90 


Evolution  Equations:  (t,  E;  Fc 


x+  =  ~(p-++p+-)x+  +P-+ 

X-  =  1  —  x+ 

Likelihoods  (Choice  1):  a,  j3  =  ± 

pap(E-,  Fc)  =  V7AGVe^)  ,  7  =  K/L’T 

T 

Activation  Energy:  cc,  (3  =  ±90,  Fc  =  2 EcPr,  AGq  = 

'  AG0(1  -  Faf}(E)/Fc )2  ,  Fa/J(£)  <  Fc 


A  GUE-FC)  = 


Taf3 

Likelihoods  (Choice  2): 

p+_(E;Ec)  = 

where 


v  0  j  Eap(E)  >  Fc 

E-+(E)  =  2  EPr  ,  F+_  =  —  F_+ 


7i 


erfcx[72(-.E  -  Ec)\ 


i  P-+(E;Ec)  = 


7i 


erfcx[72(.E  -  Ec)\ 


7i  = 


1  2Vrj 


V 


t  V  nkT 


,  72  = 


2kTrj 


Densities: 


Ka 


Vc(Fc)  —  ^  '  Otk4>k{Ec)  >  <l>k(Ec) 


=-[ln(Fc)-ln(Fc)]2/2(^)2 


o^Fc\J 2tt 


Ka 


VI {El)  =  ^  a  y  /3kVk{Ei)  ,  (fk{Ei)  =  ,  7 — e 


^^k~k0 


=  1  r-E?/2±D2 

k. 


ctJ\/2tt 


(16) 

(17) 

(18) 


(19) 


(20) 


(21) 


(22) 


(23) 


Table  2:  Components  of  the  polarization  model  based  on  180°  dipole  switching  and  the  physical 
assumptions  (4)  with  Pp>  =  P^  and  r  =  riso- 
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because  they  allow  accurate  evaluation  of  the  density  expansions  (5)  and  construction  of  lookup 
tables  to  improve  algorithm  efficiency  as  detailed  in  Section  3.2. 

If  we  let  FCi,Eij  and  Vi,Wj  respectively  denote  the  quadrature  points  and  weights,  we  obtain  the 
approximate  relations 

Ni  Nj 

xa(E  +  Ep  Fc)r'/(£,/)i/c(Fc)d£'/dFc  «  EE-  a(E  +  Ep ;  FCi )vI( Ejj )  zzc  ( Fc. ) u j Wj .  (24) 

i=l  j= 1 

For  the  validation  results  reported  in  Section  5,  sufficient  accuracy  was  achieved  using  a  trapezoid 
rule  with  Nt  =  Nj  =  41. 


3.2  Determination  of  Phase  Fractions  and  Likelihoods 

From  (24),  it  is  observed  that  evaluation  of  the  model  requires  the  solution  of  Ni  x  Nj  differential 
equations  for  each  phase  fraction  xa  which  in  turn  requires  Ni  x  Nj  evaluations  of  the  likelihoods  paf3- 
Due  to  its  relative  simplicity,  we  first  illustrate  a  highly  efficient  implementation  algorithm  for  the 
180°  polarization  model.  In  Section  3.2.2,  an  analogous,  but  slightly  more  complicated,  algorithm  is 
developed  for  the  90°  polarization-strain  model. 


3.2.1  Polarization  Model 

To  approximate  the  solution  of  (19),  it  is  shown  in  [4]  that  for  sufficiently  small  stepsizes  At,  an 
implicit  Euler  discretization  yields  the  stable  difference  relation 

r. ,fc+l  k+1  k  I  k+ 1 

a^_|_  -  O-j-  Jb _ |_  I-  L-2 

where  tk  =  kAt,  E k  =  Ee(tk )  =  E(tk)  +  Ej,  xk+  =  x+(tk]  e|,  Fc)  and 

1 


rk+ 1  _ 

Cl  — 


1  +  At\p — \-(Ek+1-,Fc)  +  P-\ — (Ek+1;  Fc)] 


rk+l  _ 

Co  — 


p.+{Ek+^Fc)At 


(25) 


1  +  A t\p.+(Ek+1-Fc)  +  p-i — (Ee+i;  Fc)] 


fc+i. 


From  (16)  and  (24),  it  follows  that 


Ek 


P(Ek)  = - PR  +  2PRVTXkW. 

V 


where  the  vectors 


V‘  = 


vivc(FCl),  •••  ,vnmc(Fc  ) 

1  J  1  x  Ni 


(26) 


(27) 


WT  =  ^wiuI(EIl),  •••  ,wNjvI(ElN.)\^N 

incorporate  the  quadrature  weights  and  densities  evaluated  at  the  quadrature  points.  The  TVj  x  Nj 
matrix  X+  quantifies  the  dipole  orientations  at  the  quadrature  points  and  is  defined  as 

*++1l  =  x+(tk+1,Ek+1+El3;FCi) 

Hj 

x+(tk,Ek  +  EL-FCi) 


1  +  At [p.+ (E*+i  +  Ei. ;  FCj )  +  p+_  (E*+i  +  ETj ;  Fc , )] 
p.+iE^  +  E^F.AAt 


(28) 


+ 


1  +  A t\p-+(EW  +  Ej-F „.)  +  p+_(£fc+1  +  FCj)} 
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Figure  3:  Structure  of  the  3-D  array  C\  and  2-D  array  C\  for  specified  field  value. 


Remark  1.  Because  the  constants  c\  and  C2  in  (25)  depend  on  field  values  but  not  explicitly  on  time, 
they  can  be  evaluated  offline  and  stored  for  highly  efficient  subsequent  model  implementation. 
To  achieve  this,  one  can  construct  3-dinrensional  arrays  C\  and  C2  whose  coordinates  are  Fc,  Ej 
and  E  with  values  at  the  quadrature  points  (FCi ,  Ey.  )  and  field  values  Ef  uniformly  distributed 
between  minimum  and  maximum  field  values  Emin,Emax  of  operation;  see  Figure  3.  During 
model  execution  —  e.g.,  during  parameter  estimation  or  model-based  control  implementation 
-  one  would  access  the  2-dimensional,  Ny  x  Nj.  arrays  C[  =  C\(:,  :,Ee)  and  Cfj  =  C-2{ :,  -,Ee), 
for  E<!-  nearest  the  specified  field  value  Ek+1,  and  update  X+  using  the  componentwise  matrix 
multiplication  and  summation 

Xl+l  =  C[.  xXk  +  Cl 

When  combined  with  (26),  the  efficiency  of  implementing  the  thermal  relaxation  model  is  the 
same  as  that  of  Algorithm  2.6.4  on  page  117  of  [18]  for  the  negligible  relaxation  model. 

An  alternative  is  to  form  arrays  V _ _ storing  the  likelihoods  at  the  Nt  x  Nj  quadra¬ 

ture  points  and  predefined  field  values.  The  matrices  C'f,  C:j  can  then  be  constructed  using 
the  operations  defined  in  (28).  This  is  slightly  less  efficient  but  can  preferable  if  memory  is 

limited  and  V _ \-,V-\ _ are  used  for  the  Jacobian  construction  described  in  Section  4.  It  is  also 

advantageous  if  variable  stepsizes  At  are  required  for  control  implementation. 


3.2.2  Polarization-Strain  Model 

From  (8)-(10),  it  follows  that  the  discretized  polarization  and  strain  relations  can  be  expressed  as 

P(E,a)  =  x°E+  Y  (ada  +  Pr)  VT Xa(t\  E,a)W 

a=±,90  ,  „ 

(29) 

e(E,  a)  =  sEa  +  Y  (Eda  +  £r)  Vt Xa(t\  E,  a)W 

a=±,9  0 

where  V  and  W  are  defined  in  (27).  The  Nt  x  Nj  matrices  Xa(E,a)  are  defined  componentwise  by 

[Xa(t;E,a)\ij  =  xa(t;E  +  EI:j,a,FCi)  ,  a  =  ±,90  (30) 

where  x±$o  are  solutions  of  (11). 

For  discrete  time,  field  and  stress  values  ty-  =  kAt,,Ek  =  E(tk)  and  ak  =  a(tk ),  we  let  = 
xa(tk ;  Ek  +  Ej,ak,  Fc)  and  [Xk\ij  =  xa(tk;  Ek  +  Ej:/ .  ak,  FCi)  denote  the  associated  dipole  fractions 
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and  matrix  of  dipole  fractions  evaluated  at  quadrature  points.  An  implicit  Euler  discretization  then 
yields 


a 


a , 


k+l  k+l 
11  x- 

k+l  k+l 
21  x- 


,k+lk+l 


*12 


X, 


X+l^k+l 


*22 


X, 


=  x_  +  a 


=  Xj 


fc+i 

13 

7fe+l 

l23 


(31) 


where 


flii"1  =  1  +  At  (p_9o  +  p — h  +  P90-)  ,  a^i1  =  ~At(P-+  +P90+) 


k+l  _ 


k+l  _ 
ll2 

k+l  _ 


=  ~At(p+-  +P90-) 


-.k+l  _ 
l22 


=  1  +  A t  (p+ go  +  p+-  +  P90+) 


<4+1  =  Atpgo+- 


*13  —  Atp90_ 

We  note  that  like  the  dipole  fractions,  these  constants  depend  on  (E,  Ej,  Fc)  ;e.g,  ak+1(E,EI,Fc)  = 
Atp90_(Ek+1  +  EriFc). 

Using  Cramer’s  rule,  one  can  solve  (31)  to  obtain 


x 


k+l 


=  C 


,fe+l 


k+ 1 


=  d 


1 

•k+l  k 


X_  +  Ex*  +  6 


Jc+1 


1 


I  +d 


k+l  k 


jk+l 


where 


rk+ 1  _ 

(  i  — 


dk+1  = 


r.k+1 

x22 

det 


ck+1  = 


—a 


k+l 

12 


det 


rk+ 1  _ 

Co  — 


—a 


fc+i 

21 


det 


4+1  = 


X+ 1 
2n 

det 


1 

det 

4+1  =  E 

3  det 


k+lk+l 
22  °13 


—  a 


fc+l „fc+l 
12  M23 


k+l  k+l 


11 


23 


—  a 


k+lk+l 
13  U21 


j  > 

min 


det  =  a\txa\tX  -  af+'a^1. 

To  facilitate  efficient  implementation,  we  follow  the  strategy  detailed  in  Remark  1  and  construct 
3-D  arrays  C\  -  V3  comprised  of  the  constants  c\  -  d,3  evaluated  at  the  quadrature  points  ( FCi ,  Ej 
and  field  values  Ef-  uniformly  distributed  between  the  minimum  and  maximum  field  values  E 
and  £max;  e.g.,  [C2]ij  =  At\p+-(Ee  +  E^-,  FCi)  +p90-{Ee  +  ETj\  FCi)]/det(£^,  EI;j,  FCi).  During  model 
implementation,  one  then  accesses  the  NixNj  matrices  C\.  C|,  63,  D[,  D l  D3,  nearest  to  the  specified 
field  value  Ek .  The  matrices  of  dipole  fractions  are  then  updated  using  the  relation 

Xk++1  =  C{.  x  Xl  +  Cl  x  Xk  +  C3 

Xk+1  =  D{.  xX^  +  De2.  xXk+De3 

\rk~\~l  T  \rk~\~l  vk-\~  1 

A90  —  1  A  +  A- 


which  again  involves  only  componentwise  matrix  multiplication  and  summation. 


4  Parameter  Estimation  Problem 

The  data  for  the  parameter  estimation  problem  is  taken  to  be  either  polarization  measurements  Pi 
or  polarization-strain  pairs  {Pi,£i}  corresponding  to  Nj  field  values  E(ti).  The  model  response  with 
parameters  p  for  the  same  field  values  and  a  constant  prestress  is  denoted  by  P(E(ti ),  ao',p)  and 
e(E(ti),a0;p). 

To  calibrate  the  180°  polarization  model,  summarized  in  Table  2,  with  the  parameter  p  given  by 
(15),  we  minimize  the  functional 

..  Nd 

J(P)  =  (32) 

i— 1 
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where 


e(U;p)  =  P(E(ti),(ro;p)  -  Pi.  (33) 

We  first  note  that  the  magnitudes  of  the  first  four  parameters  {?y,  Pr,  7,  r}  can  vary  from  1CP4  to  108 
so  the  performance  of  optimization  algorithms  is  significantly  improved  by  scaling  the  parameters 
by  dividing  by  their  magnitudes.  To  accomplish  this,  we  form  a  vector  s,  whose  components  are  the 
scale  of  each  parameter,  and  minimize  instead 


1 


Nd 


J(q)  =  ,£e2(^.  x  s) 


(34) 


i=  1 


where  q  =  p./s  and  ./  and  .x  respectively  denote  componentwise  division  and  multiplication.  Because 
the  components  of  q  are  on  the  order  of  unity,  (34)  is  more  readily  mimimized  using  standard 
optimization  software. 

To  estimate  parameters  for  the  combined  polarization  and  strain  model  (8),  we  employ  the 
functional 

1  Nd  1  Nd 

J (q)  =  7)  '^2wPep(ti-,q-  x  s)  +  -  ^2w£el(ti-,q.  x  s) 

(35) 


2  ■  '  '  •  '  '  '  2 

i=  1  i= 1 


2Nd 


where 


i=  1 

ep(U ;  q.x  s)  =  P(E(ti),a0 ;  g.xs)  -  Pi 
eE{ti ;  q.  x  s)  =  e(E(ti),  <r0;  q.  x  s)  -  e* 

_  f  wPe2P{ti\q.  x  s)  ,  *  =  1,  ■  ' 

1  w£ej(ti;q.  x  s)  ,  *  =  1X^+1 


,2JVd. 


Here  p  =  g.  x  s  is  given  by  (7).  Since  polarization  is  typically  on  the  order  of  10  1  and  strain  on  the 
order  of  10~3,  the  weights  wp  and  w£  are  chosen  to  balance  the  polarization  and  strain  components 
of  (37). 

The  behavior  of  the  functionals  J(q)  is  typical  of  those  arising  in  parameter  estimation  problems; 
they  have  shallow  slopes  near  the  global  minimum  (if  it  exists)  and  they  are  typically  nonconvex 
with  multiple  local  minima. 

One  possibility  is  to  employ  stochastic  optimization  techniques  such  as  genetic  algorithms,  sim¬ 
ulated  annealing,  and  differential  evolution  [22],  These  techniques  reduce  the  reliance  on  accurate 
initial  parameter  estimates  and,  in  theory,  provide  global  convergence.  However,  their  convergence 
rates  are  slower  —  they  may  require  infinite  time  for  convergence  —  and,  because  they  are  nonde- 
terministic,  multiple  optimizations  can  yield  varying  final  parameter  values. 

Alternatively,  one  can  employ  gradient-based  techniques  such  as  the  interior-reflective  Newton, 
Levenberg-Marquardt,  or  sequential  quadratic  programming  methods  employed  in  the  MATLAB 
routines  lsqnonlin  and  fmincon.  The  efficiency  and  success  of  gradient-based  optimization  methods 
is  predicated  on  determining  good  initial  parameter  estimates  and  being  able  to  accurately  determine 
gradients.  We  show  next  that  analytic  gradient  expressions  can  be  constructed  for  the  homogenized 
energy  model  and,  in  Section  4.2,  we  use  physical  aspects  of  the  model  to  develop  data-driven 
techniques  to  provide  initial  parameter  estimates  based  on  measured  properties  of  polarization  and 
strain  data.  In  combination,  this  permits  routines,  such  as  lsqnonlin  which  is  used  to  optimize  (32) 
and  (33),  to  be  implemented  in  a  matter  of  minutes. 
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4.1  Analytic  Jacobian  or  Gradient  Relations 

Most  commercial  optimization  packages  provide  the  option  for  constructing  the  Jacobian  using  fi¬ 
nite  difference  approximation  or  providing  an  analytic  expression.  The  former  option  is  certainly 
feasible  for  optimizing  (34)  or  (35)  but  it  will  be  significantly  slower  and  slightly  less  accurate  than 
optimization  using  the  analytic  Jacobian.  We  illustrate  first  the  computation  of  the  Jacobian  for 
the  polarization  model  (16)  with  the  functional  (34)  since  it  is  simpler  than  that  for  the  combined 
polarization-strain  model. 


Polarization  Model 

The  Jacobian  corresponding  to  (34)  is 


J  = 


de(t\\p) 

dr] 


9e(typ) 

dPR 


de(tN ,;p)  de(tN ,;p) 


dr) 


9Pr 


where  e(tp,p )  is  defined  in  (33).  We  note  that 

de(U;p)  d  d 

P(E(t.),»0;p)=- 


de(ti;p) 

90k 


de(tN.;p) 


90k 


(36) 


[J\ij  =  P(E(U),a0-,p )  =  ^-P(E(U),a  0;p)^  =  £ -P(E(U),a0-,p)Sj  (37) 

oqj 


)dp£  =  _d_ 

dqj  dpj 


so  the  Jacobian  is  simply  the  gradient  or  sensitivity  of  P  with  respect  to  the  parameters.  For  the 
discretized  model,  these  derivatives  have  the  form 


dP 

dr] 

dP 

$7 

dP 

dak 


E 

r,2 


2  PRVT[d1X+]W 


2Pr 

<*e 


(Vk-V)T[X+]W  , 


dP 

dPR 

dP 

Ih  : 

dP 

dpk 


--  -1  +  2 VT[X+]W  +  2 PRVT[dpRX+]W 

2  pRvT[dTx+\w 

VT[X+](Wk-W ) 

rLe  Pi 


where  [VT]i  =  ViUc(FCi),  [W]j  =  WjV^E^),  as  indicated  in  (27),  and  [Vp]i  =  Vi4>k(FCi),[Wk\j  = 
u>j(pk(Eij).  The  components  of  the  JVj  x  Nj  matrices  [d^X+]ij  =  ^^(i;  E(t)  +  Ekj ,  <7o;  FCi)  for 
£  =  Pr,  7,  t,  can  be  found  by  differentiating  the  evolution  equation  (19)  and  switching  the  order  of 
differentiation  to  obtain  the  differential  equations 

xPr  =  -  ( dpRp+ _  +  dpRp-+)  x+  -  ( p+ -  +  p-+)xPR  +  dpRp-+ 
x1  =  -  (<97p+_  +  <97p_+)  x+  -  ( p+ _  +  p_+)x7  +  d^p-p 


XT  =  —  {dTp+~  +  drP-  +  )  X+  ~  (p+~  +  P-  +  )xT  +  drp-  + 


in  the  variables  xPr  =  ,  x 7  =  and  xT 


dPRP+- 


-7 


(- E+Ec )2 
2EC 


P+- 


0 


E  <EC 
E  >  Ec 


dx+ 

~5E 


with 


dpRP-+ 


-7 


( e-ec )2 

2  Ec 


P-  + 


0 


E  <EC 
E>  Ec 


<97P+-  =  -A  G%_p+- 
dTp+-  =  -~p+- 

T 


,  (97p_+  =  -A  Ga_+p-+ 
,  dTp-+  =  -~p-+. 

T 


(38) 


Given  values  of  x+  and  x_,  the  differential  equations  (36)  can  be  solved  using  an  implicit  Euler 
method  analogous  to  that  described  in  Section  3.2.1  to  approximate  the  solution  of  (19). 
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Remark  2.  The  Jacobian  construction  can  be  facilitated  by  the  use  of  lookup  tables,  in  a  manner 
analogous  to  that  described  in  Section  3.2.1,  to  avoid  the  repeated  evaluation  of  the  likelihoods 

P-i _ and  p _ for  all  of  the  quadrature  points  (FCt ,  Ej. ) .  Alternatively,  if  memory  is  limited, 

one  can  construct  two  arrays  V _ \-,V-\ _ storing  the  likelihoods  at  field  values  and  quadrature 

points  in  a  manner  similar  to  that  discussed  in  Remark  1  of  Section  3.  If  memory  is  plentiful,  it 
will  be  more  efficient  to  construct  lookup  tables  associated  with  the  quantities  defined  in  (38). 


Polarization-Strain  Model 


The  construction  of  the  Jacobian  for  the  functional  (33)  is  analogous  but  slightly  more  compli¬ 
cated  since  it  involves  more  parameters  p  =  [PR,  £r,  £9r,  x°\  d+,  -sE,  7,  790,  rjgo,  oik,Pk\-  We  first  note 
that  the  model  (29)  exhibits  the  parameter  dependencies 

P(E,a)  =  XaE+  J2  C vda  +  P%)VT(ak)Xa(t-E,a,p)W((3k ) 

a=±,90 

e(E,  a)  =  sEa  +  ( Eda  +  eaR)  VT(ak)Xa(t-  E,  a,p)W(/3k ) 

a=±,90 


where  V  and  W  are  defined  in  (27),  Xa  is  defined  in  (30),  and 

P=  [P^,4.,4V+,7,^90,T180]. 

With  the  assumptions(4)  and  definitions  [Vj]«  =  Vi4>k{FCi ),  [Wk\j  =  u)j(pk(Ejj ),  it  follows  that 
BP 

=  VT\X+]W  -VT\X_]W+  T  ( ada  +  P%)VT[dp+Xa\W 

R 


bpr 

gp 

Bd+ 

BP  _ 

~BC  = 

BP 


a=±,  90 


aVT[X+}W -aVT[X_]W+  Y  K  +  PR)VT[Bd+Xa]W 


a=±,90 


Y  (<?da  +  P%)VT[dcXa\W 


a=±,90 


E 


+  PS) 


BP 


<*=±,90  ^ 

BP  „  BP 
=  E  ,  =  0 


(V-V)T[Xa]W  ,  —=  Y 


d(3k  a^±, 90  'EtP* 


(ada  +  PR)  vt^x^w  _  w ^ 


dx 


8sE 


for  (  =  £r,  e9r,  7,  T90  and  riso-  Differentiation  of  the  strain  relations  yields  similar  terms. 

We  next  discuss  the  construction  of  the  Nt  x  Nj  matrices  [8^Xa]ij  =  j^{t\E{t)  +  Ej^ao;  FCi) 
where  £  =  PR.  d+,  £R,  4°,  7, 730  and  riso-  Differentiation  of  (11)  with  respect  to  £  and  switching  the 
order  of  differentiation  yields  the  fourteen  differential  equations 


X-  =  ~^(P-90  +P-+  +P90-)X-  ~  {p-90  +P-+  +P90-)4 

d  c  dpQo- 

+  3£<J?+~  ~p"~>x+  +  w+-  ~P9 °-)x+  +  ~Qf~ 

3  (39) 

x+  =  g^(p-+  ~  P90+)x-  +  (p-+  ~  P90+)x‘L 

0  dpQo+ 

~g^(P+90  +  P+-  +  P90+)X-  -  {P+90+P+-  +P9 0+)4  +  — 7^4 
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in  the  variables  x^+  =  ,  .xl  =  For  an  input  field  value  E(tk),  x+  and  x_  can  be  determined 

using  the  techniques  detailed  in  Section  3.2.2,  and  the  solution  to  (39)  can  be  approximated  using  an 
implicit  Euler  discretization  and  Cramer’s  rule  in  a  manner  analogous  to  that  detailed  in  Section  3.2.2. 
To  evaluate  the  derivatives  of  the  likelihoods 


Pa0 


=  _jLe-7AG^  = 
Faf 3 


E—e-'yAG0(i-Fap/Fc)2  F  <  f 

Taf 3  ’  —  C 

0  ,  Fap  >  Fc 


we  note  that 


and 


dpag 

<9+ 


—  -AGapPa/3 


dpag  -1 

>  a  —  Pag 
orag  rap 


dpag 

di 


27AG0 


-FaP/Fc)paP^  ,  FaP<Fc 
0  ,  Faf 3  >  Fc. 


The 


values  of  resulting  from  the  relations  (13)  are  compiled  in  Table  3. 


90+ 

+90 

90- 

-90 

+- 

-+ 

pi 

E 

-E 

—  E 

E 

-2  E 

2  E 

4 

a 

—a 

a 

—a 

0 

0 

4° 

—a 

a 

—a 

a 

0 

0 

d+ 

Ea 

—Ea 

-Ea 

Ea 

-2  Ea 

2Ea 

Table  3:  Values  of  the  derivatives 


4.2  Data-Driven  Techniques  for  Initial  Parameter  Estimates 

An  ideal  data  set  for  obtaining  data-driven  initial  parameter  estimates  and  optimization-based  final 
parameter  values  has  two  components:  (i)  saturated  major  loop  and  biased  minor  loop  polarization 
and  strain  data,  and  (ii)  creep  polarization  or  strain  data  at  two  or  more  field  values  as  depicted  in 
Figure  4(a)  and  (b).  If  creep  data  is  not  available,  frequency-dependent  data,  such  as  that  depicted 
in  Figure  4(c)  can  be  used  to  provide  an  initial  estimate  for  rgo- 


Figure  4:  (a)  and  (b)  Saturated  major  loop,  minor  loop,  and  creep  data  used  to  determine  initial  esti¬ 
mates  for  -Pr  ,  d+,  7,  rgo  and  rigo-  (c)  Frequency-dependent  strain  data  used  to  determine 

initial  estimate  for  rgo  or  riso- 
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Figure  5:  Polarization  and  field  values  used  to  obtain  initial  estimates  for  Fc.  07,  aj. 


As  noted  in  (7),  the  parameters  to  be  estimated  are  p  =  [Pr,  ej,  e^,  xa ■>  d+,  sE,  7,  rgo,  riso,  07,  (3k\- 
The  representations  (5)  for  the  densities  vc(Fc)  and  vi(Ej)  additionally  requires  estimation  of  the 
coefficients  p  =  [Fc,  07,  07]. 

Data-driven  techniques  to  estimate  initial  values  for  ,  ej,  xa ,  d+,  Ec,  Fc,  07 ,  07  are  summa¬ 
rized  in  Table  4  with  reference  to  Figures  4  and  5.  The  values  of  07  and  07  are  then  used  to  construct 
the  vectors  07  =  {2^07}  and  07  =  {2 kac}  for  the  density  basis  functions  <pk{Ej )  and  <pk{Fc)-  One 
strategy  to  obtain  initial  values  for  the  coefficients  ag.  and  /3k  is  simply  to  take  oik  =  (3k  =  1- 
Documented  values  provide  an  initial  estimate  for  sE . 

Initial  Estimates  for  7,  rgo  and  riso  Based  on  Creep  Data 

Creep  data  collected  at  two  or  more  fixed  fields,  as  indicated  by  the  thick  lines  in  Figure  6,  can  be 
used  to  determine  initial  estimates  for  7,  rgo  and  riso-  Whereas  either  strain  or  polarization  data  can 
be  employed,  the  former  exhibits  more  sensitivity  to  90°  switching  so  we  employ  it  here.  The  fixed 
field  values  are  denoted  by  E\  and  E2  as  indicated  by  (1)  and  (2)  for  the  strains  plotted  in  Figure  6. 
We  note  that  E\  designated  by  (1)  satisfies  E\  <  —  Ec.  Finally,  let  to  designate  the  time  at  which 
the  field  was  fixed,  let  £7,62  and  £3  designate  measured  strain  values  at  equally  spaced  subsequent 
times  t\  =  to  +  At,  t2  =  to  +  2 At  and  t3  =  to  +  3 At,  and  let  e2  =  £2  —  £1,  e3  =  £3  —  £2. 


Assumption  1.  We  assume  that  the  decay  behavior  of  dipole  fractions  can  be  adequately  approx¬ 
imated  by  neglecting  interaction  fields  and  assuming  that  switches  occur  at  ±EC.  For  a  fixed 

field  E  >  Ec,  we  assume  that  X-  =  0  and  that  pgo-  =  P+90  =  P-\ _ =  0.  Similarly,  for  E  <  —Ec, 

we  assume  that  x+  =  0  and  that  p-  go  =  P90+  =  P _ b  =  0. 


With  these  simplifying  assumptions,  it  follows  from  (2)  that 

X-  =  -P90-X-  +P90-  ,  E  <  -Ec 
X+  =  P90+X+  +  P90+  ,  E  >  Ec 
For  fixed  fields,  (40)  has  the  analytic  solution 

x-{t,E)  =  1  -  (x-(t0,E)  -  i)e-(t-to)P9o-(E) 
x+(t,E)  =  1  -  (x+(t0,E)  -  l)e~ 


(40) 
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(t-to)p90+(E) 


E  <  —E, 
E  >  Ec. 


(41) 


Parameter 

Data  or  Source 

Justification 

Pr 

Remanence  polarization 

4 

Renranence  strain 

c90 

£r 

Minimum  strain 

xa 

Slope  of  E-P  curve  prior  to  switch¬ 
ing. 

From  the  relation  P  =  xaE  =  eo(er  —  1  )E, 
it  follows  that  the  relative  permittivity  er 
and  \a  satisfy  xa  =  eo(er  —  1)  where  eo  = 
8.85  x  ICR12  F/m  is  the  permittivity  of  free 
space. 

d+ 

Slope  of  E-e  curve  prior  to  switching. 

Can  be  compared  to  documented  material 
values. 

sE 

Slope  of  a-e  curve. 

Can  be  compared  to  documented  material 
values. 

7>r90>Tl80 

Creep  data  as  shown  in  Figure  4. 

7  and  rgo  given  by  solution  of  (43)  using  data 
from  two  fixed  fields  E\  and  We  typi¬ 

cally  take  tiso  =  ^r90. 

a-k  =  Pk  =  1 

A  priori  choice 

Ec 

Fc  _  = 

2PrEc 

Coercive  field  Ec 

Whereas  the  mode  Ec  most  closely  coincides 
with  Ec  (see  Appendix  B),  by  approximating 
with  the  median  Ec,  the  relation  Ec  =  e^c 
from  (50)  can  be  used  to  obtain  /ic  =  Ec. 

(Xi 

a 

Switching  begins  when  i 'i(Ej)  intersects 
vc{— Ec)  with  2cr/  specifying  the  94.5%  con¬ 
fidence  level. 

(X  c 

Ef.  Polarization  where  switching 
completes  for  increasing  E. 

Occurs  when  z//(E/)  quits  intersecting 
uc(Ec)  at  94.5%  confidence  level.  From  (51), 
this  yields  ac  =  ^  In  ((Ef  —  2crj)/Ec)  since 

eiic+2ac  +  2(JI  =  Ef. 

Table  4:  Data-driven  algorithms  for  initial  parameter  estimates  with  reference  to  Figures  4-6. 


From  Assumption  1,  it  follows  that 

s(E,  a)  «  sEa  +  E  daxa  +  £RXa 

o==b,90  a=d=,90 


which,  when  combined  with  the  physical  assumptions  in  (4),  yields 


£3 

e-2 


[( Ed+  +  —  £^°)x+(t3,  E)  +  £^°]  -  \{Ed+  +  —  £^)x+(t2,E)  +  gf;0] 

[(Ed+  +  e+  -  e^j ?)x+(t2,  E)  +  efj°]  -  [(Ed+  +  -  e9£)x+(ti,E)  +  ef°] 

X+(t3,E)  -  X+(t2,E ) 
x+(t2,E)  -  x+(ti,E) 


(42) 
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Figure  6:  Zero  prestress  strain  data  from  [24]  used  to  obtain  initial  estimates  of  7  and  tqq. 


for  E  >  Ec.  The  expression  for  E  <  —Ec  is  identical.  From  (3),  (41)  and  (42),  it  follows  that 

,  E  <  -Ec 
,  E  >  Ec. 


—e~7A  Gsa-W  = 

—  In  1 

(  e2 

T~90 

-to 

<1 

Ve3 

—e~^A  G9o+(E)  = 

—  In  1 

(  e2 

A)0 

■+0 

< 

Ve3 

(43) 


We  note  that  if  \E\  >  Ec ,  as  is  the  case  for  the  data  at  E\  indicated  in  Figure  6  by  (1),  energy 
barriers  are  eliminated  so  that  AG°Lg0(E)  =  AG°ig0(E)  =  0  which  yields 


At 


T90  = 


ln(e2/e3) 


In  this  case,  the  decay  data  at  the  second  fixed  field  can  be  used  to  solve  for  7.  If  |£,|  <  Ec  for  both 
fixed  fields,  or  if  Ec  is  unknown,  one  can  eliminate  190  to  first  solve  for  7.  To  illustrate  for  E\,  E2  <  0 
and  corresponding  strain  differences  ei2,ei3  and  e22,e23,  elimination  of  790  yields  the  relation 

e7AG-9o(®i)  ln(e12/e13)  =  e7AG“90(£’2)  ln(e22/e23)  (44) 


for  7.  Note  that  (44)  can  be  easily  solved  for  7  by  plotting  the  solution  as  a  function  of  72  and 
approximating  the  root  or  by  using  a  symbolic  routine;  e.g.,  the  symbolic  MATLAB  command 
solve. m.  The  parameter  790  is  then  determined  from  (43).  An  initial  value  of  riso  is  then  taken  to 
be  Tiso  =  t7r90- 
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The  data-driven  algorithm  to  determine  initial  parameter  estimates  for  the  180°  switching  polar¬ 
ization  model  (16)  is  analogous. 

Initial  Estimates  for  rgo  Based  on  Frequency-Dependent  Data 

Frequency-dependent  strain  or  polarization  data,  such  as  that  plotted  in  Figures  4  and  7,  can  be 
used  in  a  similar  manner  to  obtain  an  initial  estimate  for  rgo  or  nso-  We  let  to  denote  the  time  at 
which  a  field  reversal  occurs  at  Emax  or  Em ;n.  As  before,  we  let  £±,£2  and  £3  designate  measured 
strains  at  equally  spaced  subsequent  times  ti  =  to  +  At,  t2  =  to  +  2At  and  t3  =  to  +  3At,  and  let 

e2  =  £2  —  £i,  e3  =  £3  —  £2. 


Assumption  2.  For  a  field  reversal  at  Emax,  we  make  the  simplifying  assumption  that  x_  ~  0 
and  that  dipole  kinetics  are  dominated  by  switching  from  90°  to  180°  orientations;  i.e.,  pgo-  = 

p_i_go  =  p~\ _ =  0.  Similarly,  for  reversal  at  Em ;n,  we  assume  that  x+  ps  0  and  p- go  =  P90+  = 

p _ [_  =  0.  These  are  the  same  assumptions  made  in  Assumption  1.  Finally,  we  assume  that 

activation  energies  are  sufficiently  small  to  permit  the  approximation 

P90+  =P90-  =  —  ■  (45) 

790 


With  these  assumptions,  (11)  reduces  to 
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which  has  the  solution 
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x_(i,  E)  =  1  -  (x-(t0,  E)  -  l)e-(*-*°)/T90  ?  E  « 

x+(t,E)  =  l-(x+(t0,^)-l)e-('-to)/T9°  ,  E  « 


(46) 


(47) 


If  we  assume  that  E\d+  ~  E-2<i+  ~  Ej,d+  «  Em ax,  analysis  similar  to  that  used  to  obtain  (42)  yields 
the  relation 

T9°  =  ln(e2/e3)  (48) 

which  provides  an  initial,  data-driven  estimate  for  rgo-  Due  to  the  assumption  (45),  which  is  more 
stringent  than  the  case  for  creep  data  where  fields  are  fixed,  this  data  does  not  readily  admit  a 
data-driven  algorithm  to  initially  estimate  7.  Hence  it  must  simply  be  estimated  through  a  least 
squares  fit  to  data. 


Figure  7:  Frequency-dependent  strain  data  used  to  obtain  an  initial  estimate  for  rgg. 
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5  Model  Calibration  and  Validation  for  PZT  and  BaTiC>3 

We  illustrate  the  performance  of  the  data-driven  parameter  estimation  algorithms  using  PZT  data 
reported  in  [24]  and  BaTiC>3  data  from  [5]. 

For  the  PZT  data,  polarization  and  strain  measurements,  in  the  absence  of  applied  mechan¬ 
ical  stresses,  were  collected  using  a  P-802.00  piezo-actuator  whereas  a  prestressed  NEC  Model 
AE0505D08  actuator  with  a  restoring  spring  was  used  to  collect  stress-dependent  data.  Both  are 
multilayer  devices  with  similar  dimensions  so  we  briefly  summarize  the  latter  and  refer  the  reader 
to  [24]  for  details  regarding  both  devices. 

The  AE0505D08  actuator  was  comprised  of  80  active  layers,  each  of  thickness  0.1  mm.  The 
dimensions,  including  packaging,  were  6.5  x  6.5  x  10  mm  which  yielded  the  cross-sectional  area 
A  =  42.25  mm2  =  4.225  x  10~5  m2.  The  polarization  was  measured  using  a  Sawyer- Tower  circuit 
comprised  of  a  reference  capacitor  connected  in  series  with  the  PZT  actuator.  Displacements  were 
measured  using  a  fiber-optic  displacement  sensor  having  an  effective  resolution  of  0.1  /xm  and  strains 
were  computed  by  dividing  the  total  thickness  (8  mm)  of  the  active  layers. 

When  reporting  the  accuracy  of  model  fits,  we  define  the  relative  errors  to  be 

_  i  {wP[P{E(ti),cr;  q )  -  Pt\2  +  we\e{E{ti),cr\  q )  -  £;]2) 

6re*“  ESiM? +  «;.*?) 

The  functionals  (34)  and  (35),  with  the  analytic  Jacobians  constructed  using  the  algorithms  in 
Section  4.1,  were  minimized  using  the  MATLAB  nonlinear  least  squares  routine  lsqnonlin  with  the 
options  ‘Algorithm’=‘trust-region-reflective’,  ‘TolFun’=le  —  4,  ‘Jacobian' =‘on’  and  a  lower  bound  of 
0  specified  for  all  parameters  except  which  was  bounded  below  by  —1  x  10~2. 

5.1  Zero  Prestress  PZT  Data 
90°  Polarization-Strain  Model 

We  first  illustrate  the  calibration  and  performance  of  the  90°  switching  polarization-strain  model 
(8)  and  180°  switching  polarization  model  (16)  using  zero  prestress  data  reported  in  [24]  obtained 
using  a  P-802.00  PZT  actuator. 

The  data  used  to  obtain  initial  parameter  values  is  plotted  in  Figures  6  and  8  and  resulting  values, 
obtained  using  the  algorithms  in  Table  4  are  compiled  in  Table  5.  Corresponding  model  predictions 
and  initial  density  representations  are  respectively  plotted  in  Figure  9  and  Figure  11(a)  and  (c).  The 
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Figure  8:  Zero  prestress  polarization  and  strain  data  from  [24]  used  to  obtain  initial  estimates  of 
Pr  )  £R>  Xa,Fc,  crc  and  07. 
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Figure  9:  Polarization  and  strain  model  prediction  to  zero-prestress  data  from  [24]  using  initial 
parameter  values  obtained  using  the  data-driven  algorithms  in  Table  4:  (a)  time  domain  and  (b)-(d) 

phase  space;  ( - )  experimental  data,  (-•-•)  model  fit,  ( - )  model  fits  for  8  selected  loading 

cycles.  Labels  (1)  and  (2)  in  Figure  6  indicate  strain  data  used  to  find  initial  estimates  for  7  and  rgo- 


initial  values  of  e^,  and  \a  are  obtained  directly  from  the  polarization  and  stain  data  shown 
in  Figure  8.  The  density  coefficient  Fc  =  PrEc  results  from  Ec  =  0.5  x  106,  and  2<t/  =  0.5  x  106 
yields  07  =  0.25  x  106.  Similarly,  the  relations  in  Figure  5  yield  ac  =  0.51og((1.5  —  0.5)/0.5)  «  0.35. 
Finally,  the  initial  estimates  for  7  and  rgo  were  computed  using  the  data  shown  in  Figure  6.  The 
initial  value  =  374  x  10~12  is  a  literature  value  of  <^33  reported  for  PZT. 

Whereas  the  model  quantifies  the  trends  exhibited  by  the  data,  the  initial  predictions  are  not 
sufficiently  accurate  for  material  or  device  characterization.  We  then  minimized  the  functional  (35) 
using  2140  data  values  to  obtain  the  optimized  parameters  reported  in  Table  5  and  the  resulting 
fits  and  densities  shown  in  Figure  10  and  11(b)  and  (d).  The  optimized  relative  error  is  10.99%  It 
is  observed  that  the  calibrated  polarization  and  strain  models  accurately  fit  the  data  in  all  drive 
regimes  including  the  saturation  and  burst  regions,  the  creep  regions,  and  minor  loops.  Moreover, 
use  of  the  analytic  Jacobian  and  efficient  implementation  algorithms  permitted  the  optimization  to 
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-9° 

eR 

X+ 

d+ 

7 

190 

Tiso 

Init. 

0.2 

0.081x  10~2 

0.0 

4.71xl0-8 

374xl0“12 

1.38  x  10"4 

0.64 

.064 

Opt. 

0.18 

0.050X10"2 

-5.64xl0~2 

5.35X10"8 

930.46  x  10“12 

2.08xl0~3 

0.41 

.043 

Table  5:  Initial  and  optimized  parameter  values  for  the  polarization  and  strain  model  with  zero 
prestress.  The  fixed  density  parameters  were  Fc  =  0.2  x  106,(j j  =  0.25  x  106  and  ac  =  0.35. 
The  optimized  density  parameters  07  =  1.65,  ao  =  0.92,  ck_i  =  0.77,  a_2  =  2.55,  a_3  =  6.69  and 
Pi  =  0.79,  /?o  =  1 .24,  /3_i  =  1.12,  P-2  =  2.95,  P-3  =  9.45  were  obtained  using  initial  estimates  of  1. 
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Figure  10:  Optimized  fit  of  the  polarization  and  strain  model  to  zero-prestress  data  from  [24]:  (a) 

time  domain  and  (b)-(d)  phase  space;  ( - )  experimental  data,  (—•  —  •)  model  fit,  ( - )  model 

fits  for  8  selected  loading  cycles. 


be  completed  in  673  s  on  a  Mac  Pro  with  a  2.26  GHz  processor  and  8  GB  of  memory. 

We  note  that  the  final  parameter  values  in  Table  5  are  quite  close  to  the  initial  values  predicted 
by  the  data-driven  algorithms.  Hence  much  of  the  inaccuracy  in  the  initial  model  prediction  is  due 
to  the  initial  density  parameter  choices  =  1.  The  initial  values  of  and  (3k  can  be  easily 

tuned  to  provide  more  accurate  fits;  however,  the  choice  of  1  simplifies  the  data-driven  algorithms 
and  is  sufficiently  accurate  so  that  the  optimization  routine  yields  accurate  fits  in  a  timely  manner. 


x  10  ' 


(a) 


(b) 


(d) 


Figure  11:  (a)  Initial  and  (b)  optimized  critical  driving  force  density  vc{Fc)  with  41  equally  spaced 
quadrature  points  marked  as  dots,  (c)  Initial  and  (d)  optimized  interaction  field  density  vj(Ej). 


180°  Polarization  Model 

For  applications  and  operating  regimes  with  negligible  to  moderate  prestresses  <7o,  the  field- 
polarization  response  exhibits  minimal  dependency  on  90°  switching.  In  these  cases,  the  polarization 
can  be  adequately  modeled  using  the  180°  switching  model  summarized  in  Table  2. 
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Pr 

7 

790 

Init. 

0.2 

0.25x10s 

2.79xl0“3 

0.15 

Opt. 

0.20 

0.22x10s 

3.78x  10-3 

0.10 

Table  6:  Initial  and  optimized  parameter  values  for  the  180°  model  zero  prestress. 


To  illustrate,  we  employed  the  data-driven  algorithms  of  Section  4.2  to  obtain  the  parameter 
estimates  compiled  in  Table  6.  We  then  minimized  the  functional  (34)  using  2140  data  points  to 
obtain  the  optimized  fit  and  densities  shown  in  Figure  12.  The  accuracy  of  the  fit  demonstrates  the 
applicability  of  the  180°  switching  model  for  characterizing  the  polarization  in  low  to  moderate  stress 
regimes.  Moreover,  the  use  of  the  implementation  algorithms  of  Section  3  and  analytic  Jacobian 
relations  of  Section  4.1  renders  the  model  calibration  codes  highly  efficient  and  optimization  was 
achieved  in  approximately  36  s  on  the  Mac  Pro. 


(a)  (b) 

Figure  12:  Optimized  fit  of  the  180°  polarization  model  to  zero-prestress  data  from  [24]:  (a)  time 

domain  and  (b)  phase  space;  ( - )  experimental  data,  (-•-•)  model  fit,  ( - )  model  fits  for  8 

selected  loading  cycles. 


5.2  Prestressed  Actuator 

The  estimation  of  parameters  and  performance  of  the  polarization-strain  model  for  a  prestressed 
actuator  is  illustrated  in  [19]  using  data  reported  in  [24]. 

5.3  Variable  Loading  Rates 

To  illustrate  the  physical  phenomena,  and  resulting  model  response,  associated  with  variable  loading 
rates,  we  consider  data  from  [24]  that  was  collected  at  loading  rates  of  5  V/s,  50  V/s  and  500  V/s. 
Since  this  data  was  collected  using  the  P-802.00  actuator,  that  had  purely  electrical  loading,  the 
polarization  and  strains  were  modeled  using  the  relations  (8)  that  neglect  the  prestress  ao  and  spring 
and  damping  coefficients  k  and  c  associated  with  restoring  mechanisms. 

The  strain  and  polarization  data  plotted  in  Figure  13  exhibit  differing  dynamics  at  the  three 
loading  rates,  especially  in  the  burst  region  where  +180°  to  90°  and  90°  to  —180°  dipole  transitions 
are  dominant.  In  the  field-polarization  plots,  this  is  reflected  by  delayed  first-order  minor  loops  as 
indicated  by  the  curves  (l)-(3).  The  slower  90°  to  —180°  transitions,  and  resulting  lagging  effects 
at  high  loading  rates,  are  even  more  noticeable  in  the  field-strain  plots  due  to  the  degree  to  which 
stains  depend  on  90°  dipole  switching. 

To  model  this  rate-dependent  phenomena,  the  data-driven  techniques  were  used  to  compute  the 
initial  parameter  estimates  compiled  in  Table  7  along  with  the  density  coefficients  Fc  =  0.2  x  106,  ac  = 
0.347  and  07  =  0.25  x  106.  We  then  minimized  the  functional  (43),  using  the  5  V/s  and  500  V/s 
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Figure  13:  Polarization  and  strain  data  from  [24],  collected  at  (a),(d)  5  V/s,  (b),(e)  50  V/s,  and 
(c) , (f )  500  V/s. 


data,  with  wp  =  1  ,w£  =  40,  000,  to  obtain  the  optimized  parameters  in  Table  7,  the  densities  shown 
in  Figure  14,  and  the  model  fits  plotted  in  Figure  12.  We  note  that  the  model  is  predicting  the 
50  V/s  behavior. 

The  tinre-polarization  and  time-strain  plots  illustrate  that  the  model  is  quantifying  both  the 
major  and  minor  loop  behavior  at  the  three  loading  rates.  The  phase  plots  illustrate  that  while  the 
model  is  very  accurately  characterizing  the  field-polarization  behavior,  there  is  some  inaccuracy  in 
the  field-strain  fits  and  prediction.  This  may  be  due  to  limitations  in  the  relations  used  to  quantify 
90°  and  180°  switching  mechanisms.  However,  the  model  is  quantifying  the  lagged  switching  behavior 
at  500  V/s  that  is  present  in  both  the  polarization  and  strain  data.  We  note  that  improved  fits  for 
each  loading  rate  can  be  achieved  if  parameters  are  optimized  for  that  rate.  The  dynamic  nature 
of  the  model  is  illustrated  by  the  fact  that  it  adequately  characterizes  the  rate-dependent  material 
behavior  using  one  set  of  parameters.  This  makes  it  advantageous  for  material  characterization, 
device  optimization,  and  control  design  of  applications  requiring  variable  loading  rates. 


xIO5  xIO5 


Figure  14:  (a)  Critical  driving  force  density  z'c(-Fc)  and  (b)  interaction  field  density  pj(Ej)  with  41 
equally  spaced  quadrature  points  marked  as  dots. 
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Figure  15:  Optimized  fit  of  the  polarization-strain  model  to  data  from  [24]  for  loading  rates  of  (a) 

5  V/s  and  (c)  500  V/s.  (b)  Model  prediction  for  a  loading  rate  of  50  V/s:  ( - )  experimental  data, 

( - )  model  fit. 
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5.35  x  10“8 

930.46  x  10“12 

2.08  x  10-3 

0.41 
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Opt 

0.21 

0.066  x  10"2 

-0.071  x  10"2 

3.93  x  10~8 

992.26  x  10-12 

8.37  x  10“3 

1.48  x  10~9 

0.16 

Table  7:  Initial  and  optimized  parameter  values  for  the  polarization  and  strain  model  obtained 
using  loading  rates  of  5  V/s  and  500  V/s.  The  optimized  density  parameters  are  a\  =  5.02,  cko  = 
1.82  x  10-2,  a_i  =  1.24,  a_2  =  2.56,  ct_3  =  2.89  and  ft  =  4.63  x  10~2,ft  =  1.13  x  10“8ft-i  = 
4.71  x  10— 8,  /?_ 2  =  3.42  x  10-14,/3_3  =  11.01. 

5.4  Single  Crystal  BaTi03 

Here  we  illustrate  the  performance  of  the  model  for  characterizing  the  behavior  of  single  crystal 
BaTi03  using  data  from  [5].  For  the  data  considered  here,  the  applied  stress  was  parallel  to  the 
applied  field  thus  motivating  consideration  of  a  uniaxial  model.  Furthermore,  we  employ  the  grain- 
level  or  single  crystal  relations  (10). 

We  first  determined  the  initial  and  optimized  parameter  values  reported  in  Table  8  through  a  fit 
to  the  data  collected  at  prestresses  of  -0.36  MPa  and  -1.78  MPa.  The  correlation  of  the  parameter 
values  for  ,  ej,  e^0,  \+  and  cL+  to  corresponding  properties  of  the  data  can  be  directly  observed. 
We  also  note  that  rgo  >  Tigo  which  indicates  that  180°  switching  occurs  more  quickly  than  90° 
switching.  This  is  consistent  with  analogous  observations  for  PZT,  reported  in  [23],  and  explains  the 
flat  behavior  of  the  P-e  data. 

The  resulting  fit  at  -0.36  MPa  and  -1.78  MPa  and  predictions  at  prestresses  of  -0.72  MPa  and 
0  MPa  are  shown  in  Figure  16.  It  is  observed  that  the  polarization  and  strain  fits  at  -0.36  MPa 
and  -1.78  MPa  are  quite  accurate  as  are  the  polarization  predictions  at  the  other  two  prestresses. 
Whereas  the  strain  prediction  at  -0.72  MPa  is  moderately  accurate,  the  model  overpredicts  the 
strains  for  the  unstressed  crystal  where  90°  switching  is  minimal.  This  is  in  accordance  with  the 
results  reported  in  [12]  and  indicates  some  limitations  in  the  90°  switching  characterization  in  the 
absence  of  prestresses. 
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Table  8:  Initial  and  optimized  BaTi03  parameter  values  obtained  through  a  fit  to  -0.36  MPa  and 
-1.78  MPa  prestressed  data  from  [5].  The  density  representation  (14)  employed  the  optimized  value 
Fc  =  0.0145  x  106. 


6  Concluding  Remarks 

The  homogenized  energy  model  is  a  multiscale,  nricronrechanical  framework  that  quantifies  hysteresis 
and  constitutive  nonlinearities  intrinsic  to  ferroelectric  materials.  It  incorporates  mechanisms  at  mul¬ 
tiple  spatial  scales  by  combining  energy  analysis  at  the  domain  level  with  stochastic  homogenization 
techniques  to  construct  material-level  constitutive  relations.  The  multiple  timescales  associated  with 
dipole  switching  processes  are  incorporated  by  employing  Boltzmann  principles  to  quantify  transi¬ 
tion  rates  associated  with  dipole  fractions.  In  combination,  the  framework  characterizes  a  range 
of  rate-dependent  hysteretic  phenomena  associated  with  the  materials.  Furthermore,  constitutive 
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Figure  16:  Fit  of  the  polarization  and  strain  model  to  single  crystal  BaTiC>3  data  from  [5]  with  a 
prestress  of  (d)-(f)  -0.36  MPa  and  (j)-(T)  -1.78  MPa.  Model  predictions  for  prestresses  of  (a)-(c) 
0  MPa  and  (g)-(i)  -1.07  Mpa:  ( - )  experimental  data,  ( - )  model  fit. 


relations  resulting  from  the  framework  can  be  directly  employed  to  construct  distributed  models,  for 
complex  devices,  that  are  amenable  to  finite  element  implementation. 

In  this  paper,  we  presented  algorithms  that  facilitate  efficient  model  calibration  and  implemen¬ 
tation.  The  efficiency  of  the  parameter  estimation  algorithms  is  based  in  part  on  the  physical  nature 
of  the  model  which  permits  model  parameters  to  be  correlated  with  properties  of  field-polarization 
and  field-strain  data  as  shown  in  Figure  4.  Specifically,  we  presented  data-driven  algorithms  to 
determine  initial  estimates  for  P^,  s9^,  xa ,  d+,  7,  Tgo,  tiso  that  facilitate  subsequent  optimization 

through  least  squares  fits  to  the  data.  Furthermore,  this  data  can  be  used  to  determine  values  for 
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the  density  coefficients,  Ec,Fc,ac  and  07. 

Major  loop,  biased  minor  loop,  and  creep  data,  of  the  form  depicted  in  Figure  4,  represents  ideal 
data  to  be  collected  if  there  is  flexibility  regarding  the  design  of  characterization  experiments.  How¬ 
ever,  in  many  cases,  one  is  limited  to  biased  minor  loop  data  such  as  that  shown  in  Figure  1.  In  such 
cases,  one  can  employ  values  previously  reported  in  the  literature  as  initial  values  for  optimization 
routines. 

The  coefficients  1/33  and  sE  and  relative  permittivity  er,  which  can  be  used  to  compute  Xm  are 
commonly  reported  for  ferroelectric  materials  and  documented  values  can  be  employed  as  initial 
parameter  estimates.  The  remanence  polarization  Pr.  remanence  strain  £r.  and  coercive  field  Ec  are 
also  fundamental  properties  of  ferroelectric  materials  that  are  often  reported  only  in  material  char¬ 
acterization  investigations  rather  than  in  compilations  for  transducer  design.  The  model  parameters 
Pr,  £ri  )  7,  T90  and  Tigo  are  related  to  intrinsic  material  properties  that  are  critical  to  actuator 
and  sensor  design  and  characterization. 

This  motivates  the  compilation  of  material  property  and  parameter  tables  and  libraries  that 
facilitate  the  calibration  of  models  for  devices  that  employ  these  materials.  Representative  values 
of  CZ33,  sfj,  Er,  xa  from  the  literature  and  Pr,er,Ec  from  the  data  employed  here,  and  in  the  com¬ 
panion  paper  [19],  are  respectively  compiled  in  Tables  9  and  10.  Similarly,  representative  values  of 
Pt  £>R  i  7)  r90)  ri80  f°r  PZT,  PLZT  and  BaTiC>3  are  summarized  in  Table  11.  We  note  that  while 
variations  in  these  parameters  exist  for  differing  devices  and  operating  regimes,  these  representative 
values  can  be  employed  to  initialize  optimization  routines.  Furthermore,  it  is  noted  in  Appendix  A 
that  we  are  developing  a  website  of  data  libraries,  simulation  codes,  and  modeling  frameworks  for 
ferroelectric,  ferromagnetic,  and  ferroelastic  materials.  Data  properties  and  parameters  that  can  be 
used  for  material  and  device  characterization  will  be  compiled  and  updated  at  this  site. 


d33  (pC/N) 

sfi  (m2/N) 

tr  (C/Vm) 

Xa  =  offir  -  1)  (C/Vm) 

PZT-5A 

374 

1.88  x  1CT11 

1800 

1.59  x  10“8 

BaTiOs 

85.6 

8.05  x  10~12 

163 

1.47  x  10-9 

PLZT 

1188 

1.47  x  10“n 

6356 

2.99  x  10“8 

Table  9:  Material  properties  for  PZT-5A  [17],  BaTiC>3  [10]  and  PLZT  [14].  Note  that  the  free  space 
permittivity  has  the  value  eo  =  8.85  x  10~12. 


Pr  (C/m2) 

£r  (Unitless) 

Ec  (MV/m) 

PZT-5A 

0.2 

0.081  x  10-2 

0.5 

BaTiOs 

0.26 

0.67  x  10“2 

0.05 

PLZT 

0.25 

0.25  x  10-2 

0.36 

Table  10:  Remanence  polarizations,  remanence  strains,  and  coercive  fields  for  the  PZT  data  [24], 
single  crystal  BaTiC>3  data  [5]  and  PLZT  data  [14]. 


A  Available  Codes 

To  facilitate  model  validation  and  dissemination  to  the  community,  we  have  made  codes  and  data 
available  at  the  website  http://www4.ncsu.edu/~jhcrews/smart/code/pzt/index.html. 
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n 

4 

-9° 

eR 

X+ 

d+ 

7 

T90 

T180 

PZT 

0.18 

0.050  x  10"2 

-5.64  x  10”2 

5.35  x  10“8 

930.46  x  10”12 

2.08  x  10-3 

0.41 

0.043 

BaTiOa 

0.25 

0.73  x  10"2 

-0.99  x  10^2 

2.82  x  10“8 

841.89  x  10-12 

4.1  x  10~4 

0.34 

0.00041 

PLZT 

0.24 

0.13  x  10-2 

-0.03  x  10~2 

10.06  x  10"8 

1390.5  x  10-12 

9.62  x  10-4 

0.032 

0.70 

Table  11:  Representative  parameter  values  for  PZT,  BaTi03  and  PLZT.  Appropriate  units  are  given 
in  the  notation  table  in  Section  1. 


B  Lognormal  Density 

We  summarize  here  aspects  regarding  the  lognormal  density  that  are  pertinent  to  constructing  co¬ 
ercive  field  densities.  Details  can  be  found  in  [1, 11]. 

If  Y  ~  N(n,  a)  is  a  normal  random  variable  with  mean  /j  and  standard  deviation  a,  X  =  exp(Y) 
is  lognormally  distributed,  X  ~  lnN(n,cr),  with  density 

=  _J^e-(lns-M)2AA  (49) 

xayj  ^ 

The  mean  x,  mode  x,  median  x  and  variance  a2  of  X  are  related  to  [i  and  a  by  way  of  the  relations 
x  =  E[X\  =  e^2/2  ,  x  =  e^2 

/  2  \  2  (50) 

x  =  ,  u2  =  var[X]  =  (ea  -  1J  e2fl+a  . 

As  illustrated  in  Figure  17,  the  mean,  mode  and  median  differ  due  to  the  skewed  nature  of  the 
density.  We  note  that  for  the  normal  density,  n  and  a  have  the  same  units  as  x  whereas  the  physical 
meaning  of  these  units  is  lost  in  the  lognormal  relations  (50). 

The  dependence  of  /  on  [i  and  a  is  illustrated  in  Figure  18.  Note  that  X  approaches  a  normal 
density  with  x  =  x  =  x  =  efl  as  cr  decreases  and  the  density  limits  to  the  Dirac  density  at  as 

(7  — >  0. 

As  noted  on  page  9  of  [1],  the  qth  quantile  of  the  lognormal  density  can  be  expressed  as 

[e",e"]  (51) 


Figure  17:  Mode  x,  median  x  and  mean  x  of  the  lognormal  density  with  n  =  0,  a  =  1. 
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as  illustrated  in  Figure  19.  This  relation,  along  with  (50),  will  be  used  in  Section  4.1  to  specify  initial 
density  parameter  values  in  terms  of  measured  properties  of  polarization  and  magnetization  data. 
For  example,  two  standard  deviations,  which  represents  a  probability  of  94.5%  is  plotted  for  normal 
and  lognormal  densities  in  Figure  19. 


(a)  (b) 

Figure  19:  Median  values  and  94.5%  confidence  intervals  with  n  =  1,  a  =  0.5  for  a  (a)  normal  density 
[n  —  2u,  /i  +  2<r],  and  (b)  lognormal  density  [eM_2<T,  eM+2<7] . 
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